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Toward  immersed  boundary  simulation  of  high 
Reynolds  number  flows 

By  Georgi  Kalitzin  and  Gianluca  Iaccarino 


1.  Motivation  and  Background 

In  the  immersed  boundary  (IB)  method,  the  surface  of  an  object  is  reconstructed  with 
forcing  terms  in  the  underlying  flow  field  equations.  The  surface  may  split  a compu- 
tational cell  removing  the  constraint  of  the  near  wall  gridlines  to  be  aligned  with  the 
surface.  This  feature  greatly  simplifies  the  grid  generation  process  which  is  cumbersome 
and  expensive  in  particular  for  structured  grids  and  complex  geometries. 

The  IB  method  is  ideally  suited  for  Cartesian  flow  solvers.  The  flow  equations  written  in 
Cartesian  coordinates  appear  in  a very  simple  form  and  several  numerical  algorithms  can 
be  used  for  an  efficient  solution  of  the  equations.  In  addition,  the  accuracy  of  numerical 
algorithms  is  dependent  on  the  underlying  grid  and  it  usually  deteriorates  when  the  grid 
deviates  from  a Cartesian  mesh. 

The  challenge  for  the  IB  method  lies  in  the  representation  of  the  wall  boundaries  and 
in  providing  an  adequate  near  wall  flow  field  resolution.  The  issue  of  enforcing  no-slip 
boundary  conditions  at  the  immersed  surface  has  been  addressed  by  several  authors  by 
imposing  a local  reconstruction  of  the  solution.  Initial  work  by  Verzicco  et  al.  (2000) 
was  based  on  a simple  linear,  one-dimensional  operator  and  this  approach  proved  to  be 
accurate  for  boundaries  largely  aligned  with  the  grid  lines.  Majumdar  et  al.  (2001)  used 
various  multidimensional  and  high  order  polynomial  interpolations  schemes.  These  high 
order  schemes,  however,  are  keen  to  introduce  wiggles  and  spurious  extrema.  Iaccarino 
& Verzicco  (2003)  and  Kalitzin  & Iaccarino  (2002)  proposed  a tri-linear  reconstruction 
for  the  velocity  components  and  the  turbulent  scalars.  A modified  implementation  that 
has  proven  to  be  more  robust  is  reported  in  this  paper. 

The  issue  of  adequate  near  wall  resolution  in  a Cartesian  framework  can  initially  be 
addressed  by  using  a non-uniform  mesh  which  is  stretched  near  the  surface.  In  this  paper, 
we  investigate  an  unstructured  approach  for  local  grid  refinement  that  utilizes  Cartesian 
mesh  features. 

The  computation  of  high  Reynolds  number  wall  bounded  flows  is  particularly  chal- 
lenging as  it  requires  the  consideration  of  thin  turbulent  boundary  layers,  i.e.  near  wall 
regions  with  large  gradients  of  the  flow  field  variables.  For  such  flows,  the  representation 
of  the  wall  boundary  has  a large  impact  on  the  accuracy  of  the  computation.  It  is  also 
critical  for  the  robustness  and  convergence  of  the  flow  solver. 


2.  Numerical  Technique 

The  computational  code  IBRANS  is  based  on  the  steady-state  incompressible  RANS 
equations  closed  either  with  an  one-equation  turbulence  model  (Spalart  & Allmaras 
1994)  or  a two-equation  model  (Kalitzin  & Iaccarino  2002).  A second-order,  implicit, 
cell  centered  Cartesian  discretization  is  used  within  a SIMPLE  pressure- velocity  coupling 
algorithm  with  a segregated  solution  of  the  field  equations  (Ferziger  & Peric  2002). 
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The  equations  are  discretized  in  three-dimensions  with  a typical  7-point  stencil.  In 
this  paper,  however,  we  will  refer  for  simplicity  to  the  corresponding  two-dimensional 
equations.  For  each  cell  the  discretized  field  equation  for  momentum,  turbulent 

scalar  and  pressure  correction  can  be  written  as: 

+ aB<pE+1  + «+1  + a!v^jV+1  + a^S!  = S?i,j)  C2-1) 

where  the  indices  E,  W,  S,  N denominate  the  neighboring  cells  that  have  a common  face 
with  cell  (i,j),  e.g.  E for  cell  (i  + 1,  j)  etc.  The  superscript  n indicates  the  iteration.  The 
equation  is  under-relaxed  implicitly  to  increase  the  diagonal  dominance  of  the  implicit 
matrix.  The  Strong  Implicit  Procedure  (SIP)  by  Stone  (1968)  is  used  for  the  solution  of 
(2.1). 

The  flow  solver  is  set  up  as  a virtual  wind  tunnel.  The  ^-direction  is  the  main  flow 
direction  with  inflow  and  outflow,  the  other  directions  have  either  the  usual  wall  bound- 
ary, symmetry  or  periodicity  conditions.  The  domain  is  partitioned  in  x-direction  and 
the  equations  are  solved  in  parallel  using  MPI.  An  object  is  placed  into  the  flow  domain 
and  the  immersed  boundaries,  described  in  the  next  section,  are  applied  to  the  velocity 
components  and  the  turbulence  variables.  There  is  no  special  treatment  of  the  pressure 
at  the  immersed  boundary.  The  immersed  boundaries  affect  the  Poisson  equation  for  the 
pressure  correction  only  through  the  source  term  which  represents  the  divergence  of  the 
velocity  field.  Thus,  the  Poisson  equation  is  solved  in  the  entire  computational  domain. 
A parallel  geometric  multigrid  is  being  used  to  accelerate  the  convergence  of  the  Poisson 
equation. 

The  non-uniform  computational  grid  and  the  IB  interpolation  stencil  is  generated 
automatically  and  the  pre-processing  time  to  start  the  computation  is  negligible.  The 
procedure  is  based  on  the  ray  tracing  technique  that  allows  one  to  identify  the  cells  cut 
by  the  immersed  boundary.  An  initial  coarse  (and  typically)  uniform  mesh  is  specified 
by  the  user;  a three-steps  iterative  procedure  is  employed:  (a)  tag  the  grid  identifying 
the  cells  cut  by  the  immersed  boundary,  (b)  split  these  cells  in  the  direction  that  moves 
the  cell  center  closer  to  the  immersed  boundary,  (c)  regenerate  a structured  grid  by 
propagating  the  cell  split  to  the  domain  boundaries,  and  restart  from  (a).  This  procedure 
is  repeated  until  a prescribed  distance  from  the  wall  is  achieved  for  all  the  interface  cells. 
The  geometry  definition  is  based  on  the  StereoLitography  format  and  therefore  the  CAD 
model  can  be  used  directly. 


3.  Immersed  Boundary  Treatment 

After  the  pre-processor  detects  the  computational  cells  that  are  cut  by  the  body,  the 
cells  are  divided  in  those  that  are  inside  and  outside  of  the  body  as  shown  in  Fig.  1. 
The  cut  cells  are  separated  in  two  type  of  cells  corresponding  to  the  location  of  their 
cell  center.  The  cells  that  have  their  cell  center  outside  or  on  the  surface  are  labeled  as 
interface  cells.  The  other  cut  cells  are  treated  like  inside  cells. 

The  flow  variable  <j>,  which  represents  a velocity  component  or  a turbulence  variable,  is 
set  to  zero  in  the  inside  cells.  Note  that  the  turbulence  models  considered  have  zero  wall 
boundary  conditions  for  all  turbulence  variables.  At  the  interface  cells,  the  nearby  wall 
is  modeled  with  an  off  wall  boundary  condition  which  in  general  consists  of  an  implicit 
linear  interpolation  stencil  with  an  explicit  (non-linear)  correction. 

A sketch  of  the  linear  interpolation  algorithm  is  shown  in  Fig.  la.  Like  the  plotted  tan- 
gential velocity  component,  the  scalar  variable  cj>  in  the  interface  cell  (i,  j)  is  interpolated 
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Figure  1.  Schematic  representation  of  interpolation  algorithm  for  the  tangential  velocity.  Gen- 
eral case  (left)  and  wall  function  approach  for  a grid  aligned  case  (right).  Tangential  velocity: 
— > computed,  — > interpolated  linearly  on  intermediate  location,  — > imposed  in  interface  cell. 


using  only  the  neighboring  cells  that  are  further  away  from  the  body.  The  interpolation 
is  carried  out  in  two  steps.  First,  for  each  neighboring  cell  an  interpolation  is  carried  out 
normal  to  the  wall  to  an  intermediate  location  that  has  the  same  distance  to  the  wall  as 
the  interface  cell  center  This  interpolation  might  be  linear,  quadratical  or  other, 

depending  on  the  near  wall  behavior  of  the  considered  variable.  The  second  interpolation 
is  on  a surface  that  is  parallel  to  the  wall.  Here  we  use  the  inverse-distance  weighted 
method  proposed  by  Franke  (1982)  that  has  the  property  of  preserving  local  maxima 
and  producing  smooth  reconstructions.  The  interpolation  of  in  the  interface  cell  is: 


<%,.,•)  = YlWm4>m/q  with Wm  = ’ q = Wm  and  p = 2 t3-1) 


where  <t>m  represents  the  previously  interpolated  values,  hm  is  the  distance  between  the 
location  of  4>{i,j)  and  the  location  of  4>m  and  H represents  the  maximum  of  the  hm's. 
The  sum  is  over  all  neighboring  cells  with  m = W,  E,  S,  N.  The  weights  wm  are  zero 
for  all  non-qualifying  neighbors.  For  each  interface  cell  the  described  linear  interpolation 
results  in  a set  of  weights  /3m  which  correspond  to  the  effect  of  the  neighboring  cells  on 
the  interface  cell: 

- Pw<l>w  + Pe<I>e  + + Pn4>n  (3-2) 

By  considering  only  the  neighboring  cells  ( W, , E,  S,  N)  that  have  a common  face  with 
the  interface  cell  (i,j)  the  linear  interpolation  (3.2)  can  easily  be  treated  implicitly; 
the  implicit  matrix  of  the  discretized  equation  (2.1)  can  be  modified  with  the  weights 
fjm  without  introducing  new  diagonals  with  non-zero  elements.  The  inclusion  of  other 
neighbors  such  as  corner  cells  in  the  interpolation  stencil  (3.2)  can  be  done  explicitly. 
For  a non-linear  near  wall  behavior  of  <f>  an  explicit  correction  is  added  to  the 

right  hand  side  of  equation  (3.2).  This  correction  is  computed  as  the  difference  of  the 
non-linear  value  of  and  the  linear  interpolated  value,  both  computed  with  values 
from  the  current  iteration. 

The  described  interpolation  can  be  applied  to  the  velocity  vector  v without  coupling 
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implicitly  the  momentum  equations  for  the  Cartesian  velocity  components.  For  this,  the 
velocity  components  are  transformed  into  a wall  aligned  2-D  coordinate  system  with 
the  axis  pointed  in  the  tangential  t and  normal  n velocity  direction:  vn  — (v  ■ n)n  and 
vt  = v — vn,  respectively.  The  normal  unit  vector  is  obtained  from  the  wall  distance  d as: 
n = Vd/|W|.  Both  components  vt  and  vn  are  interpolated  according  to  their  physical 
behavior  normal  to  the  wall  and  using  the  interpolation  (3.1)  along  the  surface.  The 
tangential  and  normal  velocity  components  in  the  interface  cell  are  then  rotated  back 
into  the  base  Cartesian  system.  Because  of  linearity  the  weights  (5  remain  the  same  for  the 
Cartesian  components.  As  before,  the  non-linearity  is  treated  as  an  explicit  correction. 

For  low  Reynolds  numbers,  the  near  wall  behavior  of  the  tangential  velocity  is  linear. 
As  plotted  in  Fig.  la,  the  tangential  component  on  the  intermediate  location  for  cell  N is: 
{vt)N  — ( vt)Nd(ij)/d yv-  The  corresponding  quadratically  interpolated  normal  component 
is:  (vn)N  = (wn) For  high  Reynolds  numbers,  the  near  wall  behavior  of  the 
tangential  velocity  is  in  general  non-linear  and  the  value  at  the  intermediate  location  is 
determined  with  a non-linear  procedure  that  utilizes  adaptive  wall  functions. 

In  Kalitzin  et  al.  (2003),  an  adaptive  wall  function  concept  has  been  developed  for 
body  fitted  grids  and  zero  pressure  gradient  flow.  These  wall  functions  are  independent 
of  the  location  of  the  first  grid  point  above  the  wall.  A look-up  table  provides  an  explicit 
dependency  of  the  friction  velocity  on  the  local  Reynolds  number.  This  look-up  table 
consists  of  cubic  splines  that  approximate  piecewise  a solution  obtained  on  a very  fine 
grid  for  a zero  pressure  gradient  boundary  layer  flow.  This  universal  law  is  written  in 
terms  of  the  local  Reynolds  number  Rey  rather  than  in  terms  of  y+ . The  local  Reynolds 
number  is  defined  as  Rey  = U+y+  = Uy/v.  The  look-up  table  is  turbulence  model  specific 
as  the  universal  law  varies  slightly  depending  on  the  turbulence  model.  Similar  look-up 
tables  provide  the  explicit  dependency  of  the  eddy-viscosity  and  turbulence  variables 
on  the  local  Reynolds  number  or  on  the  wall  distance  y+.  The  tables  also  include  the 
derivative  of  the  velocity  and  turbulence  variables. 

In  this  paper,  we  only  consider  a high  Reynolds  number  case  in  which  the  plate  is 
aligned  but  not  coincident  with  the  grid  lines  as  shown  in  Fig.  lb.  For  zero  pressure 
gradient  flows,  the  Navier-Stokes  equations  simplify  near  the  wall  to: 

, ,dU  o 

{}'  + ''')lhj  =Ut  (3>3) 

In  the  current  approach,  the  velocity  in  cell  (i,j)  is  chosen  such  that  the  flux  at  the  face 
( i,j  + |)  fulfills  equation  (3.3).  The  tangential  velocity  and  the  wall  distance  in  cell  N 
define  the  local  Reynolds  number  Rey  which  is  related  to  the  friction  velocity  uT  via  the 
adaptive  wall  function  table.  The  friction  velocity  uT  and  the  distance  between  the  face 
(i,j  + |)  and  the  wall  allow  to  introduce  a wall  distance  y+  which  determines  an  eddy- 
viscosity  value  The  discretized  equation  (3.3)  defines  the  tangential  velocity 

in  (i,j): 

v? 

u(i,j)  =uN-  ■ T (dN  - d(iJ))  (3.4) 

v + W(ij+i) 

This  approach  can  not  simply  be  generalized  for  a situation  shown  in  Fig.  la  as  the 
tangential  velocity  in  (i,  j)  is  not  uniquely  determined  by  one  single  interface  flux.  For  the 
general  non-aligned  case  we  employ  therefore  a crude  approach  by  simply  extrapolating 
of  the  tangential  velocity  component  and  eddy-viscosity  using  wall  function  information 
about  the  derivative  in  cell  N. 
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Figure  2.  Velocity  magnitude  for  immersed  boundary  layers  simulations  on  non-aligned  grids, 
Re  = 1000,  mesh  with  every  8th  line  plotted. 


Immersed  boundary  results  for  a flat  plate  boundary  layer  at  Re  — 1 , 000  are  shown 
in  Fig.  2.  The  three  simulations  are  carried  out  on  uniform,  non-aligned  grids  with  a 
different  angle  between  the  plate  and  the  gridlines.  The  flow  is  laminar  and  y+  of  the 
first  cell  center  is  about  0.1.  The  case  of  10°  is  the  most  challenging  as  the  wall  distance  of 
the  first  grid  point  above  the  wall  fluctuates  by  moving  in  the  streamwise  direction  along 
the  plate.  There  are  about  30  cells  in  the  boundary  layer  and  20  cells  across  the  plate 
resolving  the  cylindrical  leading  edge.  The  solution  is  nearly  independent  of  the  angle 
between  plate  and  gridlines  as  shown  with  the  velocity  magnitude  contours  in  Fig.  2d 
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and  Fig.  2e.  The  contour  lines  for  larger  velocity  values  deviate  slightly  for  the  45°  what 
is  likely  to  be  a numerical  dissipation  effect.  However,  the  skin  friction  along  the  plate 
matches  remarkably  well  the  body  fitted  solution  shown  in  Fig.  3a. 

In  Fig.  3b,  skin  friction  is  reported  for  an  immersed  boundary  simulation  of  a flat  plate 
boundary  layer  at  Re  — 10®.  The  plate  is  aligned  with  the  grid  and  the  turbulent  flow 
is  computed  with  the  Spalart-Allmaras  turbulence  model.  The  wall  distance  y+  of  the 
first  cell  center  above  the  wall  is  about  25.  The  results  obtained  using  wall  functions  are 
compared  to  the  ones  obtained  using  the  simple  linear  interpolation  and  to  a body  fitted 
solution.  The  skin  friction  levels  for  the  linear  interpolation  are  incorrect  whereas  the 
wall  functions  results  are  in  good  agreement  with  the  body  fitted  solution  downstream 
of  £ = 0.25.  Upstream  of  this  location,  the  discrepancy  between  the  wall  function  and 
body  fitted  result  is  related  to  the  grid  resolution  and  wall  function  implementation.  The 
resolution  of  the  cylindrical  leading  edge  is  very  coarse  for  this  Reynolds  number.  The 
grid  is  the  same  as  the  one  used  for  the  simulations  in  Fig.  2. 


4.  Local  grid  refinement 

A new  version  of  IBRANS  with  a local  grid  refinement  capability  is  in  development 
for  an  efficient  clustering  of  cells  in  the  boundary  layers.  The  present  implementation 
is  an  extension  of  the  “classical”  adaptive  mesh  refinement  (AMR)  technique  to  allow 
non-isotropic  refinement.  It  can  also  be  interpreted  as  a generalization  of  the  procedure 
used  for  building  coarse  grids  for  geometric  multigrid  on  structured  meshes. 

The  basic  idea  was  introduced  in  Durbin  & Iaccarino  (2002)  for  a finite  difference 
discretization.  The  AMR  grid  is  considered  as  a coarsened  version  of  an  underlying, 
structured  grid;  on  this  underlying  grid  the  cells  are  defined  (in  two-dimensions)  by  a 
couple  of  vertices  with  indices  (i,j)  and  (i  + l,j  + 1).  On  the  AMR  grid  each  element 


Immersed  boundary  simulations  375 

is  bounded  by  the  gridlines  passing  through  the  vertices  (i,j)  and  ( i + A i,j  + Aj). 
The  effective  element  size  in  the  AMR  grid  is  not  constant  (Aj  ^ Aj  and  both  indices 
depend  on  ( i,j )).  Therefore,  the  cells  are  not  organized  in  a structured  way  with  one-to- 
one  neighbors  in  each  Cartesian  direction.  This  requires  a modification  of  the  algorithm 
to  deal  with  hanging  nodes.  In  Durbin  & Iaccarino  (2002)  this  was  simply  based  on  a 
second-order  interpolation.  In  the  current  finite  volume  implementation,  flux  conservation 
is  enforced  and  the  algorithm  resembles  the  unstructured  face-based  algorithm  described 
in  Ferziger  & Peric  (2002)  and,  more  closely,  the  AMR  discretization  used  in  Ham  et 
al.  (2002).  Each  face-flux  is  computed  using  the  two  adjacent  cells  and  for  each  cell 
the  fluxes  (in  general  more  than  four)  are  collected  to  build  the  corresponding  diffusive 
and  convective  operators.  In  two  dimensions,  the  implicit  discretization  yields  a sparse 
matrix  with  elements  not  organized  in  five  diagonals  as  for  its  structured  counterpart. 
This  complexity  in  the  matrix  structure  prevents  the  use  of  the  SIP  procedure  and  a 
Krilov-type  algorithm  with  a simple  Jacobi  pre-conditioner  was  implemented.  Standard 
conjugate  gradient  is  used  for  the  pressure  equation  and  the  BiCGStab  (Van  den  Vorst 
1992)  for  the  momentum  and  turbulent  scalars. 

The  major  advantage  of  the  present  approach  with  respect  to  classical  OCTREE- 
based  (Berger  & Aftosmis  1998)  and  fully-unstructured  (Ham  et  al.  2002)  schemes 
lies  in  the  economy  and  flexibility  of  storing  and  retrieving  connectivity  information  due 
to  the  underlying  grid.  In  particular,  only  N cells  are  effectively  defined  on  a Ni  x Nj 
underlying  grid  and  they  are  defined  by  the  two  couples  (i,j)  and  (i  + A,,  j + Aj).  The 
total  storage  cost  is  AN  integers.  In  addition,  an  array  of  integers,  is  defined 

on  the  fine  grid  to  store  the  correspondence  between  the  underlying  cell  and  the  actual 
AMR  element.  In  other  words,  all  the  underlying  cells  included  in  the  range  (i,  i + A,-  — 1) 
and  (j,  j + Aj  - 1)  are  tagged  using  the  AMR  cell  number.  The  total  storage  required  is, 
therefore,  Nt  x Nj.  The  connectivity  information  for  each  cell  are  retrieved  consistently  to 
a structured  framework  by  indirectly  querying  the  array  ID^j) . The  neighbors  of  a AMR 
cell  are  and  7-D(i+A.+ljfc)  for  k ranging  between  j and  j + Aj  — 1,  in  the  positive 

and  negative  i-direction  respectively.  It  is  evident  that  the  approach  handles  multiple 
hanging  nodes  for  each  cell  and,  eventually,  allows  to  reconstruct  additional  connectivity 
information  without  any  increase  in  storage;  for  example  it  is  straightforward  to  identify 
all  the  vertex-based  neighbors. 

The  generation  of  AMR  grids  is  carried  out  by  creating  the  underlying  (fine)  grid  as 
discussed  before  and  then  coarsening  it  in  the  region  away  from  the  immersed  boundary. 
The  advantage  of  this  approach  is  that  all  the  cell  tagging  (ray  tracing)  can  be  performed 
on  a structured  grid  taking  full  advantage  of  the  alignment  of  the  cell  centers  and  the 
grid  nodes.  The  coarsening  and  the  generation  of  the  connectivity  information  is  the  last 
step  of  the  grid  generation  process. 

An  example  of  the  application  of  this  procedure  is  shown  in  Fig.  4 where  both  the 
underlying  and  the  AMR  grids  are  reported;  note  that  the  AMR  grid  contains  only  9% 
of  the  underlying  cells.  A sample  calculation  has  been  performed  on  both  grids  and  it  is 
reported  in  Fig.  5 in  terms  of  turbulent  kinetic  energy.  The  solution  on  the  AMR  grid 
is  remarkably  smooth  and  consistent  with  the  one  obtained  on  the  structured  mesh.  It 
converges  significantly  faster  as  reported  in  Fig.  6.  This  is  mainly  due  to  a reduced  aspect 
ratio  of  the  cells  away  from  the  immersed  boundary.  In  terms  of  computational  cost  the 
savings  is  about  70%. 


Figure  4.  Computational  grid  for  the  flow  around  an  airfoil.  Structured,  underlying  grid  on 
the  left  with  every  second  grid  line  plotted  and  locally  refined  grid  on  the  right. 


Figure  7.  Turbine  blade  with  internal  cooling  passages  (a)  and  planar  cuts  through  locally 
refined  grid:  (b)  cut  near  the  tip,  (c)  cut  near  mid-span 

The  local  grid  refinement  algorithm  has  been  implemented  in  three-dimensions  and 
it  works  well  even  with  very  complex  geometries  as  the  real  turbine  blade  shown  in 
Fig.  7.  The  blade  is  twisted  and  the  internal  cooling  system  of  the  blade  consist  of 
extremely  intricate  passages  with  multiple  serpentines,  ribs,  tip  ejection  holes,  film  cooling 
holes,  leading  edge  impingement,  etc.  Planar  cuts  through  the  locally  refined  grid  are 
shown  in  the  same  figure.  A comparison  between  both  cuts  reveal  a remarkable  consistent 
refinement  in  all  internal  flow  regions.  The  local  grid  refinement  also  follows  the  twisted 
outer  surface  of  the  blade. 


5.  Conclusions  and  future  work 

The  paper  presents  details  of  a Cartesian  immersed  boundary  method  for  RANS  flow 
simulations.  It  focuses  on  the  IB  implementation  for  mean  flow  and  turbulence  variables. 
The  immersed  boundary  is  represented  with  an  implicit,  linear,  off  wall  interpolation 
with  an  explicit,  non-linear  correction.  An  extension  of  the  method  with  adaptive  wall 
functions  is  presented  for  grid  aligned  cases.  The  paper  also  addresses  the  issue  of  flow 
resolution  at  the  immersed  boundary  and  discusses  a local  grid  refinement  algorithm. 

Skin  friction  distributions  are  presented  for  a boundary  layer  over  a flat  plate  that 
is  non-aligned  with  the  gridlines  at  Re  = 1,000.  The  results  are  nearly  independent  of 
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the  angle  between  the  plate  and  the  gridlines.  For  the  grid  aligned  case,  skin  friction  is 
presented  for  a flat  plate  boundary  layer  at  Re  — 106. 

Future  work  consists  in  continuing  the  development  of  IBRANS  and  combining  the 
described  technologies  into  one  mature  computational  code.  The  code  is  in  the  process  of 
being  extended  for  conjugate  heat  transfer  computations.  The  energy  transport  equation 
has  been  implemented  and  it  needs  the  development  of  appropriate  heat  transfer  im- 
mersed boundary  conditions.  The  adaptive  wall  functions  for  immersed  boundaries  need 
to  be  extended  for  general  grids  and  heat  transfer  problems.  It  is  planned  to  parallelize 
and  adapt  the  multigrid  procedure  for  the  flow  solver  with  local  grid  refinement. 
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